RNA_velocity

Part 4

RNA velocity

library(velocyto.R)
library(SeuratWrappers)
library(Seurat)
source("tianfengRwrappers.R")

preprocess <- function(loom_path, n){
  velodata <- ReadVelocity(file = loom_path)
  
  #匹配两次的barcode
  tfunc <- function(s)
  {
    s <- strsplit(s,".*:",fixed = F)[[1]][2]
    s <- strsplit(s,"x",fixed = T)[[1]]
    s <- paste(s,as.character(n),sep = "-")
    return(s)
  }
  
  velodata[["spliced"]]@Dimnames[[2]] = as.character(lapply(velodata[["spliced"]]@Dimnames[[2]],tfunc))
  velodata[["unspliced"]]@Dimnames[[2]] = as.character(lapply(velodata[["unspliced"]]@Dimnames[[2]],tfunc))
  velodata[["ambiguous"]]@Dimnames[[2]] = as.character(lapply(velodata[["ambiguous"]]@Dimnames[[2]],tfunc))
  return(velodata)
}

ctrl <- readRDS("ctrl.rds")
ctrl_velo <- as.Seurat(x = preprocess("./HSD0524.loom",1))
ctrl_velo <- subset(ctrl_velo, cells = WhichCells(ctrl))

ctrl_velo <- ctrl_velo %>% 
    PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
    SCTransform(vars.to.regress = "percent.mt", verbose = F,assay = "spliced") %>% 
    RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% 
    FindClusters(resolution = 0.1)


# for(i in levels(Idents(ctrl))){
#   (DimPlot(ctrl_velo, reduction = "umap",cells.highlight = WhichCells(ctrl, expression = `celltype` == i), pt.size = 0.5, label = T)+ggtitle(i)) %>% print()
# }

ctrl_velo@reductions[["umap"]] <- ctrl@reductions[["umap"]]
ctrl_velo@reductions[["umap"]]@cell.embeddings <-
  ctrl@reductions[["umap"]]@cell.embeddings[colnames(ctrl_velo),]

DimPlot(ctrl_velo, reduction = "umap",group.by = 'seurat_clusters', pt.size = 0.5, label = T)

ctrl_velo$Classification1 <- ctrl$celltype
Idents(ctrl_velo) <- ctrl_velo$Classification1
ctrl_velo@assays[["RNA"]] <- ctrl@assays[["SCT"]]

# saveRDS(ctrl_velo,"ctrl_velo.rds") 
# ctrl_velo$Classification1 <- Idents(ctrl)

ident.colors <- colors_list[1:length(x = levels(x = ctrl_velo))]
names(x = ident.colors) <- levels(x = ctrl_velo)
cell.colors <- ident.colors[Idents(object = ctrl_velo)]

names(x = cell.colors) <- colnames(x = ctrl_velo)

png("velocity.png")
show.velocity.on.embedding.cor(emb = Embeddings(object = ctrl_velo, reduction = "umap"), vel = Tool(object = ctrl_velo, slot = "RunVelocity"), n = 200, scale = "sqrt", cell.colors = ac(x = cell.colors, alpha = 0.5), 
    cex = 0.8, arrow.scale = 3, show.grid.flow = TRUE, min.grid.cell.mass = 0.5, grid.n = 40, arrow.lwd = 1, 
    do.par = FALSE, cell.border.alpha = 0.1)
dev.off()

h5ad_output <- "hsd_velo.h5ad"  

library(reticulate)
library(Matrix)
writeMM(t(ctrl_velo@assays$RNA@counts), file='combined.mtx')
writeMM(t(ctrl_velo@assays$spliced@counts), file='spliced.mtx')
writeMM(t(ctrl_velo@assays$unspliced@counts), file='unspliced.mtx')
write.csv(rownames(ctrl_velo@assays$spliced@counts), file = "genes.csv", row.names = FALSE)
write.csv(ctrl_velo@reductions$umap@cell.embeddings, file = "umap.csv", row.names = FALSE)
write.csv(ctrl_velo@reductions$pca@cell.embeddings, file = "pca.csv", row.names = FALSE)

write.csv(colnames(ctrl_velo@assays$spliced@counts), file = "cells.csv", row.names = FALSE)
write.csv(ctrl_velo@meta.data, file = "meta.csv", row.names = FALSE)
source_python('~/scRNAseq/vascular-analysis/build.py')
build(h5ad_output, pca = TRUE, umap = TRUE)

filename = "hsd_velo.h5ad"  
import anndata
import pandas as pd
from itertools import chain
#import scvelo as scv

ad = anndata.read_mtx('spliced.mtx')
ad.layers['spliced'] = anndata.read_mtx('spliced.mtx').X
ad.layers['unspliced'] = anndata.read_mtx('unspliced.mtx').X
ad.obs = pd.read_csv('meta.csv',header=0)
ad.obsm['X_umap'] = pd.read_csv('umap.csv',header=0).values
ad.obsm['X_pca'] = pd.read_csv('pca.csv',header=0).values


ad.var_names = pd.Index(list(chain.from_iterable(pd.read_csv('genes.csv', header=0,dtype='string').values.tolist()))).astype(str)
ad.obs_names = list(chain.from_iterable(pd.read_csv('cells.csv', header=0).values.tolist()))

ad.write(filename)

file.remove('combined.mtx')

[1] TRUE

file.remove('spliced.mtx')

[1] TRUE

file.remove('unspliced.mtx')

[1] TRUE

file.remove('genes.csv')

[1] TRUE

file.remove('cells.csv')

[1] TRUE

file.remove('umap.csv')

[1] TRUE

file.remove('pca.csv')

[1] TRUE

file.remove('meta.csv')

[1] TRUE

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.